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O. 

CN . Stochastic models of diffusion with excluded-volume effects are used to model many 

§"| biological and physical systems at a discrete level. The average properties of the 



population may be described by a continuum model based on partial differential 
equations. In this paper we consider multiple interacting subpopulations/species and 
study how the inter-species competition emerges at the population level. Each indi- 
vidual is described as a finite-size hard core interacting particle undergoing Brownian 
motion. The link between the discrete stochastic equations of motion and the con- 

^ ■ tinuum model is considered systematically using the method of matched asymptotic 

O. 

expansions. The system for two species leads to a nonlinear cross-diffusion system for 
each subpopulation, which captures the enhancement of the effective diffusion rate 
due to excluded-volume interactions between particles of the same species, and the 
diminishment due to particles of the other species. This model can explain two alter- 
native notions of the diffusion coefficient that are often confounded, namely collective 
diffusion and self-diffusion. Simulations of the discrete system show good agreement 



Q\ ' with the analytic results. 
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I. INTRODUCTION 



Stochastic models describing how interacting individuals give rise to collective behavior 
have become a widely used tool across disciplines — ranging from biology to physics to social 
sciences.-"- Despite their conceptual simplicity, particle-based models can be computation- 
ally intractable for large systems of interacting particles, as is often the case in practical 
applications. In such continuum population-level description based on partial dif- 

ferential equations that can capture the overall population density becomes attractive. The 
challenge is then to predict the correct population-level description of the key attributes at 
the particle level (such as interactions between individuals and evolution rules). 

In particular, the model of diffusive particles with hard-core repulsive (or steric) interac- 
tions is relevant to many applications, such as colloidal systems,- ion transport,- 1 ^ diffusion 
through polymers,- biological cell populations^ and animal behavior— In our previous 
work 11 we considered the diffusion under an external force of N hard spheres in a bounded 
domain Q e M. d , d = 2, 3 (of typical nondimensional diameter one). The particles were taken 
to be identical with nondimensional diameter e < 1 and diffusivity D. Starting from the 
particle-level description, we used a method based on matched asymptotic expansions for a 
small but finite volume fraction to obtain the continuum model. The result is a nonlinear 
diffusion equation for the one-particle probability density function p(x, t), 

^,t) = V x -{DV x [p+(N-l)ae d p 2 ]-{(x)p} in Q, (1) 

where a = 2(d — l)ir/d and f(x) is the nondimensionalized force (drift). Since the nonlinear 
term is positive for N > 1, we find that excluded- volume effects enhance the overall collective 
diffusion rate. 

The case of several different types of particle is relevant in many practical problems but 
has been paid far less attention by the mathematical community so far.— For example, 
multiple populations of interacting agents appear in traffic flow with heterogeneous agents,- 
pedestrian or animal motion (e.g. ants going in opposite directions^) and cellular tumor 
invasion (cancer and normal cells).— They are also important in ion transport through mem- 
brane channels, as in many applications ions can be heterogeneous.— Another application 
is the extreme case in which one of the populations is motionless and blocks the motion of 
the others; for example, anomalous diffusion in cell membranes due to obstruction (from, 
e.g, the membrane skeleton mesh, fixed proteins, or lipid rafts).— >^ 
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Much of the effort to describe multiple species with size-exclusion processes has been 
directed at on-lattice models, in which the motion of particles is restricted to taking place 
on a lattice and one defines certain hopping rules between lattice sites to account for the 
particle motion and interactions. A common approach is to take a continuum limit of the 
discrete model and obtain a partial differential equation (PDE) describing the average oc- 
cupancy of the agent population.— For example, this strategy applied to a multi-species 
motility model based on simple exclusion process with drift (or bias) leads to a system of 
nonlinear advection-diffusion equations.— 1 ^ More complicated rules have also been consid- 
ered; namely, particles that can bind to sites and interact not only with each other but also 
with a confined channel-like domain,— or myopic agents (in which the hopping probability 
depends on the number of unoccupied nearest-neighbor sites) .— Recently, Penington et. al. 
have considered a generalization of these specific on-lattice models to incorporate general 
interactions, and derived the associated continuum models systematically.— 

Here, we are interested in off-lattice models, where particles each undergo a continuous 
Brownian motion. We extend our previous work for identical hard-core interacting particles 
to the case when two types of particles are present. We call these two species the blues and 
reds after Ref. |l2|. Specifically, we allow for each subpopulation to have a different number 
of particles with different sizes and diffusivities and to be under a different external force. 

This modeling of such system of interacting particles is typically based on a microscopic 
approach using iV-coupled Brownian motions, where iV is the total number of particles, or 
a macroscopic approach using partial differential equations (PDEs). For two species, the 
latter consists of a system of PDEs for the two subpopulation densities. As mentioned 
before, microscopic models generally require many computationally intensive simulations 
to gain understanding of population-level behavior, and can become impractical to use. 
This is why continuum models are a very useful tool, and there is a lot of interest in 
predicting the correct macroscopic description of the particle-level attributes. This can be 
a very challenging task, especially when nontrivial interparticle interactions are present in 
the system. This is why continuum models are often defined phenomenologically (that is, 
written directly at the continuum level rather than derived from their discrete counterparts) 
or by making assumptions that cannot be related to individual behavior. For instance, 
closure approximations (which assume independence between individuals at some stage) are 
commonly used, yet often generate errors in the resulting continuum model. In this work we 



use instead a systematic approach based on the method of matched asymptotic expansions 
to derive the macroscopic model of the two species system, generalizing the result (JT|). 

The first part of this work is concerned with the derivation of the continuum PDE model 
from the microscopic particle-based model. We first introduce the microscopic description 
of the system based on ci/V-coupled stochastic differential equations (SDEs), where d = 2, 3 
is the problem dimension, and its associated Fokker-Planck (FP) equation for the joint 
probability density of the system. We then perform a systematic asymptotic analysis of 
the FP equation in the limit of small but finite particle volume fraction, which results in 
a nonlinear cross- diffusion system of PDEs for the two subpopulations densities of each 
species. The nonlinear terms ctnsG £is db result of the excluded-volume effects in the system. 
We compare solutions to this model with stochastic simulations of the microscopic model 
to assess the validity of our approach, and find very good agreement. We also compare 
solutions for finite-size particles with those corresponding to point particles to investigate 
the importance of excluded-volume effects. 

Second, we examine how the continuum model can be used to determine the transport 
coefficients. The fact that the model keeps two distinct densities enables us to distinguish 
between two alternative notions of diffusion coefficient: the collective diffusion coefficient, 
which describes the evolution of the total density, and the self-diffusion coefficient, which 
describes the evolution of a single tagged particle.-^ This is in contrast with the one-species 
model,— from which we could only "extract" the collective diffusion coefficient. The reason 
for that is that, although both diffusion coefficients are defined for a system with identical 
particles, to describe the self- diffusion coefficient we need to tag one individual particle and 
keep its probability density "accessible" in the continuum model. With the two species 
model this can easily be done by taking the tagged particle to constitute the second species. 
To our knowledge this is the first continuum PDE model that can be used to explain both 
types of diffusion. This makes it well suited to interpret experimental data from diffusion 
measurement experiments, which can often produce unexpected/misinterpreted results.-^ 

In the third part of this work, we explore the properties of the cross-diffusion PDE model. 
We find that rewriting the system in terms of its free energy and the mobility matrix can 
be a very useful tool to study the equilibria and stability of the system. In this alternative 
form, know as the gradient flow form, the evolution of the system can be interpreted as a 
probability flow down the gradients of the free energy. While this gradient flow structure is 
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relatively well understood in the scalar case (one species), the task of obtaining an associated 
free energy can be very challenging in the case of systems.— For instance, in our case we are 
only able to write down an explicit free-energy functional under some conditions. Another 
use of the gradient-flow structure is that it gives an explicit upper bound on the particle 
volume fraction for the validity of our asymptotic model. 

II. TWO SPECIES MODEL 

Our starting point is a system of N hard spheres (or disks) diffusing and interacting in 
a bounded domain Q in ~R d of typical dimensionless volume of order one. Suppose there are 
Nj) blue particles of diameter e& and constant diffusion coefficient D b and N r red particles of 
diameter e r and constant molecular diffusion coefficient D r , with iVj, + N r = N. Note that 
we could have chosen to nondimensionalize so that one of the two dimensionless diffusion 
coefficients D b or D r is set to one. However, we choose deliberately not to do so so that the 
resulting model is symmetric upon exchange of the blue and red labels. Also note that we 
are not making explicit any relationship between the molecular diffusivity and the particle 
size (as it might exist if, say, the Stokes-Einstein relation holds in the system). We assume 
that the particles occupy a small volume fraction, so that N b e b + N r e d <C 1. We suppose that 
the only interaction between particles is hard core repulsion (so that the particles cannot 
overlap), neglecting any electrostatic or hydrodynamic interaction forces. 

We denote the centers of the particles by Xj(i) G Q at time t > 0, where 1 < i < N. 
Each centre evolves according to the overdamped Langevin SDE 



where the Wj are N independent d-dimensional standard Brownian motions and f b and f r 
are the external forces on the blue and red particles, respectively. We restrict ourselves to 
the case where the external force acting on the ith particle depends only on its position 
Xj, that is, ffc = ffc(Xj) : Q — > M. d . This excludes external forces such as electromagnetic 
and hydrodynamic forces, in which case f; would depend on the positions of all the particles 
X = (Xi, . . . , Xjy). We suppose that the initial positions Xj(0) are random and, within the 
same species, identically distributed. 



dXi(t) = f 6 (Xi(t)) dt + y/2D b dWi(t), l<i<N b , 



(2a) 
(2b) 



dXi(t) = f r (Xi{t)) dt + y/2D^dWi(t), N b + l<i<N, 
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Let P(x, t) be the joint probability density function of the N particles. Then, by the Ito 
formula, P(x, t) evolves according to the linear Fokker-Planck partial differential equation 
(PDE) 

dP 

¥ = Vr [DV,P-F(x)P], (3) 

where and V^- respectively stand for the gradient and divergence operators with re- 
spect to the iV-particle position vector x G Q N . Here D = diag(D b , . . . , D b , D r , . . . D r ) 
is the diffusivity matrix and F(x) = (fb(xi), . . . , ffe(x7vj, ^(x^+i), . . . , f r (xjv)) is the dN- 
dimensional drift. Splitting the position vector x into the position vector for the blue 
particles, x b = (x 1; . . . , XjvJ, and for the red particles, x r = (x^+i, • • • , Xjy), equation (j3J) 
can be rewritten as 

dP _ 

— = V Sb ■ [D b V Sb P - F b (x b ) P] + W Sr ■ [D r V Sr P - F r (x r ) P] in fif , (4a) 

— # —* 

where F b (x b ) is the drift acting on the blue particles [first dN b components of F(x)}, and 
analogously for F r (x r ). Because of excluded- volume effects, ( l4"a|) is not defined in fl N but 
in its "hollow form" = Q N \ B e , where B e is the set of all illegal configurations (with at 
least one overlap), 

B e = {x e n N : 3i ^ j s.t. - x, || < i(e< + e,-) } , 

where e» = e b for i < N b and e» = e r otherwise. The domain of definition is known 
as the configuration space. On the collision surfaces dQ^ we have the reflecting boundary 
condition 

[V 2 P - F(x) P]-n = on dOf , (4b) 
where n G S^ -1 denotes the unit outward normal. The initial condition on P is 

P(x,0) =P (x). (4c) 

Since all the particles within the same species (blue or red) are identically distributed, the 
initial distribution Pq(x) is invariant to permutations of the particles labels within the same 
species. The form of (J4|) then means that P itself is invariant to permutations of the blue 
or red particle labels for all time. 

our goal is to reduce the high- dimensional PDE model (111) to a low dimen- 



As in Ref. 
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sional PDE model for the marginal density function of one particle. However, as mentioned 



in the introduction, instead of obtaining an equation for p(x±,t) = f P(x,t) dx 2 ---dxAr 
comprising the collective effect of iV identical particles, now we will have two marginal den- 
sity functions, one representative of the blue particles and one representative of the red 
particles. Because all blue particles are identical and all red particles are identical, we are 
interested in the marginal density function of, say, the first blue particle and the last red 
particle, given by 

6(x, t) = / P{x, t) 5(x — Xi) dx, r(x, t) = / P(x, t) <5(x — xjv) dx, 

respectively. We aim to reduce the high-dimensional PDE for P (j4"l) to a low-dimensional 
system of PDEs for b and r through a systematic asymptotic expansion as e&, e r — > 0. 



A. Point particles 

In the particular case of point particles (e& = e r = 0) the model reduction is straightfor- 
ward. In this case the N particles are independent and the domain is = Q N (no holes), 
which implies that the internal boundary conditions in (14b j) vanish. Therefore 



N b N 

P(f,t) = Jj6(xi,t) J] r(x,,t), (6) 

i=l i=N b +l 

and the evolution equations for the one-particle density functions b and r follow from inte- 
grating equation (jlall multiplied by 5(x — xi) and <5(x — xjv), respectively, over the configu- 
ration space Q N using (jSJ) 

^(x, t) = V x • [D b V x b - f 6 (x) b] , (7a) 
^(x,t) = V x ■ [D r V x r - f r (x)r] , (7b) 



in f2 x M + . The boundary conditions (14b j) become simple no-flux boundary conditions on 
the domain walls, 

0=[D 6 V Xl 6-f 6 (x 1 )6].n 1 , (7c) 
0=[D r V xl r-^.(x 1 )r].fii, (7d) 

on dfl x M + , where hi is the outward unit normal to dfl. The system is supplemented by 
initial values 

6(x,0)=6o(x), r(x,0) = r (x), (7e) 

in fl. Here &o( x ) = Jqn Pq{x, t) 5(x — xi) da?, and similarly for r . 



B. Finite-size particles 



When €b and/or e r are greater than zero, the internal boundary conditions in (14b j) mean 
the particles are no longer independent. The integration of OH) over x 2 , . . . , x^v results in 
integrals over the collision surfaces, on which P must be evaluated. When the particle 
volume fraction is small, the dominant contributions to these collision integrals correspond 
to two-particle collisions. This implies that if the case of N = 2 can be solved it can 
easily extended to arbitrary N. u For two species three types of interaction with N = 2 are 
possible: either two blue particles, either two red particles, or one blue particle and one red 
particle interacting. We note that the first two types involving two identical hard spheres 



have already been computed in Ref. 
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yielding equation ([T]). Therefore, here we only need 
to consider the third case, that is, Nf, = N r = 1. 

For one blue particle at position x x and one red particle at position x r , equation ( I4"a|) 
reads 

OP 

— (xi, x 2 , t) = V X1 • [ZW X1 P - f b (xi) P] + V X2 • [D r V X2 P ~ fr(x 2 ) P] , (8a) 

for (xi,x 2 ) G f^, and the boundary condition (14b j) reads 

[AV X1 P - f 6 ( Xl ) P] ■ ni + [A.V X2 P - f r (x 2 ) P] ■ n 2 = 0, (8b) 

on Xj G dfl and ||xi — x 2 || = ^(es + e r)- Here = / 1 1 1 1 , where is the component of 
the normal vector n corresponding to the ith particle, n = (n 1; n 2 ). We note that Hi = 
on x 2 G dQ, and that fii = — n 2 on ||xi — x 2 || = e. It is convenient to introduce e& r as the 
distance at contact between one blue particle and one red particle, €b r = (e& + e r )/2. 

We proceed to obtain an equation for 6(x, t) from (JHJ). We denote by f2(xi) the region 
available to particle 2 (the red particle) when particle 1 (the blue particle) is at xi, namely, 
fi(xi) = f2 \ B ebr (xx). Since the domain dimensions are much larger than the particle 
diameters, the volume |0(xi)| is constant to leading order. Integrating Eq. ( |8al) over fi(xi) 
yields 

db 

— (x 1? t)= V xl .[D 6 V xl 6-t(xi)6] 
at 



+ / [f 6 ( Xl ) P - 2D b V Xl P - D b V X2 P] ■ n 2 dS 2 (Q) 

+ / [AVx 2 P-fr(x 2 )P]-n 2 d5 2 . 



8 



The first integral in (j9]) comes from switching the order of integration with respect to x 2 
and differentiation with respect to xi using the transport theorem; the second comes from 
using the divergence theorem on the derivatives in x 2 . Using (18bj) and rearranging we find 

or n 

— (xi, t) = V X1 • [D fe V Xl b - f b (xi) b] - 2D b / V Xl P ■ n 2 dS 2 

at JdB tbr ( Xl ) 

+ [ {{D r - A,) V Xl P + P - f,(x 2 )]} • n 2 dS 2 (10) 

J9B ebr ( Xl ) 

= V X1 • [D b V Xl b - f 6 ( Xl )6] - D 6 / (V X1 P + V X2 P) • n 2 d5 2 . 

We denote the integral above by 

X 6r ( Xl ) = -D b [ ( V X1 P + V X2 P) ■ n 2 d5 2 . (11) 

At this stage, it is common to use a closure approximation.—^ For example, the classical 
closure approximation is to assume that particles are not correlated, that is, P(x 1; x 2 ,t) = 
6(xi, t)r(x 2 , t). However, the pairwise particle interaction — and therefore the correlation 
between their positions — is exactly localized near the collision surface dB ebr (pci). Here we 
use a systematic alternative method based matched asymptotic expansions 26 to compute 
X 6r (xi). 

C. Matched asymptotic expansions of P 

We suppose that when two particles are far apart (||xi— x 2 || ^> e br ) they are independent, 
whereas when they are close to each other (||xi— x 2 || ~ e br ) they are correlated. We designate 
these two regions of configuration space the outer region and inner region, respectively. 

In the outer region we define P oltt (xi, x 2 , t) = P(xi,x 2 ,t). By independence, we have 
that 

P out (xi, x 2 , t) = g fe (xi, t)g r (x 2 , t), (12) 

for some functions (fe(x, t) and <3V(x, t). It is important to note that these functions will 
not be the same in general since P is not invariant to a switch of blue and red particle 
labels. Also note that we could introduce more terms in the asymptotic expansion for P ou t- 
However, the subsequent analysis shows that the value of the integral X bT is invariant to the 
first-order correction of P out which we thus not need to consider further. 
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In the inner region, we set Xi = xj and X2 = xi+eb r x and define P(x 1; x, t) = P(xi, x 2 , t). 
Inserting this change of variables into (JHJ) yields 

BP 

e 2 br — (± 1 ,±,t) = (D b + D r )VlP 

+ e 6r V* ■ { [f 6 ( Xl ) - f r ( Xl + e 6r x)] P - 2D h V^P) ( 13a ) 
+ e^AVt 1 P-V* 1 -[f 6 (xOP]}, 



with 



* ' V ^ = n T~n * ' { A^P + + e 6r x) - f 6 ( Xl )] P, } (13b) 

L>b + JJ r 



on ||x|| = 1. Note that now (113bj) contains the no-flux boundary condition on the contact 



between the two particles, but not on dQ. As pointed out above, this is because we can 
assume that Xi is not close to dQ; the boundary only affects higher-order terms. In addition 
to ffT3bl the inner solution must match with the outer solution P out as ||x|| — > 00. Expanding 
( TT2l) in inner variables gives 

P(xi,x,t) -> g b (xi,t)g r (xi + e 6r x) 

(13c) 

~ g 6 (xi, t)g r (xi, t) + e 6r g 6 (xi) x ■ V il g r (x 1 ) H , 

as ||x|| — » 00. We now look for a solution to ( {TBI) of the form P(xi,x, t) ~ P^°)(xi,x, t) + 
e^P^^Xi, x, t) + • • • . The leading-order inner problem is simply 

v?p(°) = 0, 

x-VxP (0) =0 on ||x|| = 1, (14) 

P^ ~ ?b(xi, t)g r (xi, t) as ||x|| — >■ 00, 

which is trivially solved by 

P (0) (x!,x,t) = gb(xi,t)g P (xi,t). (15) 
At C(eb r ) (TT3T) reads, using ([15"]) and Taylor-expanding f b and f r , 

V?P (1) = 0, 

x-VxP«=x-A(x 1; t), on ||x|| = l, (16) 

P (1) ~ x ■ B(xi,t), as ||x|| ->■ 00, 
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with 



A(5M) 



1 



{D b V kl {q b q r ) + [f r (x!) - f b (xi)] q b q r }, 



D b + D r 



(17) 



B(xi,t) = q b V^q r . 



The solution to f|T6l) is 



p(i)( : 



xi, x, t) = x • A(xi, t) + x + 



(d-l)||x 



x 



) 



•[B(x 1 ,t)-A(x 1 ,t)]. 



(18) 



Combining (|15p and (|18[) . the inner solution is 



(19) 



x x ■ <^ [ffc(xi) - f r (xi)]g fe g r + D^V^qr - D b q r V ±1 q b \ + 0{e 2 br ) 



where q b = q b (5ti,t) and q r = g r (x 1 ,t). Comparing the inner solution ( !T9|) in the case of 
two distinguishable particles with that of two identical particles, for which P(x l5 x, t) ~ 
q 2 + egx • V^g, we find that two extra terms contribute to the inner solution now: one is 
due to the difference in the drift force acting on each particle, [f&(xi) — f r (xi)] q b q r , and the 
other owing to the different initial conditions and/or diffusivities, D^V^r — D b q r V '^q b . 
Note however that, even in the case of (physically) identical particles, we could use the two- 
species distinction to have the particles subdivided into two groups, each group with different 
initial conditions. In that case, all the physical parameters 6j, Di and fj are identical, and 
the distinction between the two groups in the model is contained in the outer functions q b 
and q r . 

Using the inner solution (fT9~|) . we can now evaluate the integral Ibr(xi) in ( TTTj) . Expressing 
it in terms of the inner variables, we obtain 



Now we use the normalization condition on P to determine the outer functions q b and q r . 
We find that q b (x) = 6(x) + 0(e br ) and g r (x) = r(x) + 0(ef r ), which will allow us to express 
Z&r(xi) in terms of the densities b and r. 



Xbr(xi) ~ e d hT 




(20) 
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D. System of equations for b and r 



Inserting the expression for Z^xi) obtained in (|20|) into ({TO]) we find, to 0(ef r ). 

— ( Xl ,t) = V Xl - (AV Xl 6-f 6 ( Xl )6 

+ etAj^feV^r - 7&rV Xl 6} + e d brlb [f 6 ( Xl ) - f r (xi)] fir), 

where 

2tt [(d - 1) A + dA] 2tt A 

Pb = -} Fi , Fi > 7fe 



(21) 



d D b + D r ' ' d D b + D r 
Equation (12T]) gives the evolution of 6 for a system with one blue particle and one red 

particle. The extension from one to N r red particles is straightforward up to 0(e br ), since 

at this order only pairwise interactions need to be considered. For N r arbitrary, the blue 

particle has N r blue-red inner regions, one with each of the N r red particles; hence there 

are N r copies of the nonlinear term in ( 12 ip . Similarly, for N b arbitrary the blue particle can 

have blue-blue pairwise interactions with any of the N b — 1 remaining blue particles; the 

contribution of a pairwise interaction between identical particles is found in ([TJ) . Thus the 

blue marginal density function satisfies 

— (x, t) = V x • (AVx& - f 6 (x)6 + (N b - l)e d b D b abV x b 

dt V n (22a) 

+ N r e d br {D b (f3 b bV x r - lb rV x b) + lb [f 6 (x) - f r (x)] br) 



in Q x ]R + . A similar procedure shows that the red marginal density r satisfies 
— (x,t) = V x • ( £> r V x r - f r (x)r + (N r - l)ef.D r arV x r 

(22b) 



+ N b ei{D r (/3 r rV x 6 - lr bV x r) + lr [f r (x) - f 6 (x)] 6r}) , 

in Q x M + . This system is complemented with no-flux boundary conditions on dQ x M + and 
initial conditions 

6(x,0) = 6o(x), r(x,0)=r (x) (22c) 

in Q. The coefficients in the equations (i = b and j — r and vice versa) are as follows 

_ 2(d-l)7r _ 2 1 [{d-l)D 1± dD J \ _ 2tt A r _ 9#n 

d ' A " d A + Dj ' 7i ~ d A + Dj ' l ^ dj 

for d = 2 or 3. 

We have obtained a nonlinear cross- diffusion system for the blue and red marginal prob- 
ability densities, which captures the enhancement (diminishment) of the effective diffusion 
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rate, due to excluded- volume interactions between particles of the same species (of the other 
species). Namely, the diffusion of one blue particle is enhanced by collisions with the other 
blues, and reduced by collisions with red particles [terms with +a6V x 6 and — 7f,rV x 5, re- 
spectively, in f)22ap ]. This will allows us to distinguish between two alternative notions of 
diffusion coefficient: the collective (or mutual) diffusion coefficient and the self-diffusion co- 
efficient (see Sec. IIIII) . Note that the reduced model (1221) has now a nonlinear drift term. 
This is effectively a "drag" term due to the different drift velocities of the red and blue 
particles. 

It is reassuring to find that we can recover the system for a single species from the system 
for two species (I2"2"j) . When the two species are identical, that is, D b = D r , e b = e r , and 
fb = f r , the governing equations (I22a|) and (I22bl) of the densities b and r are the same. Then, 
if the initial densities f)22cl) are equal, &o( x ) = r o( x ) : = Po( x ), then b(x, t) = r(x, t) := p(x, t) 
for all times. Consequently, we can replace b and r by p in fl22al) and recover the one-species 



equation ([I]), by noting that, when D b = D r , f3 b — 7 6 is equal to a [see f!22dj) ]. In this 
scenario, the nonlinear terms in ([22aj) rearrange to 

(N b - l)e d D b a P V xP + N r e d D b ((3 b - lb )pV^p = (N b + N r - l)e d D b apV x p, 

which coincides with the nonlinear term in the one-species equation (Tj[|). 

Finally, note that we have not specified any relation between the diffusion coefficients of 
a single particle (D b and D r ) and the particles' diameters (e b and e r ), even though a relation 
may exist between these parameters. For instance, if we are modeling a system for which 
the Stokes-Einstein relation holds, then the diffusivity should be inversely proportional to 
the particle's diameter. Thus it may be that not all of the four parameters D b , D r , e b , and 
e r can be chosen independently. 



E. System in matrix form 



The system (T2"2"l) may be rewritten in the form 



DlliW. | '] -F(6,r) f 
r \r 



(23a) 
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where D(fe, r) = j s the diffusion matrix, 

, (D b [l + (N b -l)e d b ab-N r e d brlb r] D b N r e d br (3 b b 
D(6,r) = | , (23b) 

\ D r N b ei(3 r r D r [l + (JV r -l)e?ar - iV^Tr^ 

and F(6, r) is the drift matrix 

, , W iV r et7 6 [fr(x) - f 6 (x)] A 

F(6,r) =| , (23c) 

N b et lr [f b (x) - f r (x)]r f r (x) / 

with the coefficients given in ( I22dj) . In order to focus on the diffusion part of ( 123a|) . we set 
F = for the rest of this section. 

An important consideration is that the reduced continuum model we have obtained for 
the collective or population-level behavior is different to the corresponding continuum limit 
of the discrete on-lattice counterpart model-^ the coefficients of our advection-diffusion 
system (123]) derived from the off-lattice model do not agree with those derived from the 



on-lattice model [cf. Eq. (3.1)-(3.2) in Ref. yjj. In particular, the lattice model does not 



include the negative terms in the diagonal entries of D (which are important to explain 
self- diffusion as we shall see below) and the advection term F is linear and does not include 
the difference f b — f r . It would be interesting to study whether hopping rules can be given 
on a discrete lattice which produce the model fl23|) at a continuum level. 

F. Numerical simulations 

The particle-based description (TSJ) of the problem, consisting of dN— coupled stochastic 
differential equations (SDE), is used as a benchmark solution to test the validity of the re- 
duced continuum model ( 122]) for the marginal densities b and r. To this end, we compare the 
solution of ( 122]) (obtained with the method of lines using a second-order finite-difference dis- 
cretization for the spatial derivatives, in the spirit of the positivity-preserving finite-difference 



scheme proposed in Ref. 1271 ). with Monte Carlo (MC) simulations of the 2N— coupled SDE 



2J) in two dimensions (d = 2). For the MC simulations, the reflective boundary conditions 



on dfl are implemented as in Ref. |28|, namely, the distance that a particle has travelled 
outside the domain is reflected back into the domain. The particle-particle overlaps are 
implemented similarly as follows. The difference e— ||Xj(t + At) — X 3 -(t + At)|| corresponds 
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to the distance that particles have penetrated each other illegally. Then we suppose that 
each particle has travelled the same illegal distance, and we separate them accordingly along 
the line joining the two particles' centers. To test the importance of steric interactions, we 
also compare with the corresponding solutions with e b = e r = 0. 

Figure [T] shows the results of a time-dependent simulation with f&(x) = f r (x) = 0, 
tt = [-1/2, 1/2] 2 , N b = 100, N r = 300, t b = 0.01, e r = 0.02, and D b = D r = 1. Initially, 
the blue particles are uniformly distributed, &o( x ) — 1, and the red particles are normally 
distributed in the x\-axis with zero mean and standard deviation 0.1, ro(x) = f(x; 0, 0.1 2 ), 
where f(x;fi,a 2 ) is the one-dimensional Gaussian (truncated and normalized so that its 
integral over Q is one). The figures correspond to time t = 0.05 and the simulation time 
step is set to At = 10 -5 . 




FIG. 1. Marginal densities b(x,t) ([I-IV]-b) and r(x,t) ([I-IV]-r) at time t = 0.05 with initial 
data 6o( x ) = 1 uniform and ro(x) normally distributed in the x-axis, f;, = f r = 0, D b = D r = 1, 
N b = 100 and N r = 300. (I) Solutions 6(x, t) and r(x, t) of Eq. ([7]) for point particles (e& = e r = 0). 
(II) Histograms for e b = e r = 0. (Ill) Solutions 6(x, t) and r(x, t) of Eq. ([22]) for finite-size particles 
{e b = 0.01, e r = 0.02). (IV) Histograms for e& = 0.01, e r = 0.02. Histograms computed from 10 4 
realizations of Eq. ([2]) with At = 10~ 5 . All four plots on the left and right have respectively the 
same color bar. 

The theoretical predictions for both point and finite-size particles compare very well with 
their simulation counterparts, while steric effects are clearly appreciable even though the 
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volume fraction of particles is only 0.102. The initial uniform density of the blue particles 
does not change in time when size-exclusion effects are ignored [Fig. [1] (I-b)], while it does 
diffuse towards the domain edges x ± 1/2 when they are not [Fig. [T] (III-b)] due to the 
non-uniform density of red particles. More details on this effect will be given in Sec. IIHI 
On the other hand, the red particles' initial profile, in which particles are concentrated in 
the central band, spreads faster when excluded-volume effects are included [Fig. [1] (III-r)] 
than when they are not [Fig. [1] (I-r)], indicating that the overall collective diffusion of the 
red species is enhanced. 



III. DIFFUSION COEFFICIENTS IN A TWO-COMPONENT SYSTEM 

There exist several characterizations of diffusion in a system of finite-size particles. In 
the case of a single species, there is the collective diffusion coefficient, which describes the 
evolution of the total concentration of the species, and the self-diffusion coefficient, which 
describes the evolution of a single tagged particle.— For two or more species, a third coeffi- 
cient, the cross-diffusion, expresses the motion of one species under a concentration gradient 
of the other species. 29 



A. Collective diffusion 



The collective diffusion coefficients (also known as main or principal diffusion coeffi- 
cients' 30 ) are the diagonal entries of the diffusion matrix D in f l23bl) . For instance, the 
collective diffusion of the blue species is 

D bb {b{^ t), r(x, t)) = D b [l + e d b {N b - l)a&(x, t) - t d br N rlb r{^ t)] . (24) 

The first term is the diffusion of a free blue particle. The second term, which is proportional 
to the excluded-volume created by the blue particles, always enhances the collective diffusion. 
This enhancement is due to biasing the random walk — in a gradient of b you are more likely 
to move towards the low-density region. In contrast, the third term (related to the excluded- 
volume created by the red particles) reduces the collective diffusion of a blue particle. As 



expected, setting N r = in IH)) yields the col 
of finite-size particles, see Eq. (13) of Ref. 
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ective diffusion coefficient for a single species 
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Now, with the two species in play the 



collective diffusion coefficient displays a compromise between the enhancement due to the 
finite-size interactions within your own species and the diminishment due to the presence of 
the other species. 

B. Cross-diffusion 

The cross- diffusion coefficients are the off-diagonal entries of the diffusion matrix D in 
(I23bl) . which are always non- negative. For instance, the cross-diffusion coefficient of a blue 
particle across the red particles is 

D br (b(x,t)) = D b e d br N r f3 b b(^t). (25) 

This term represents a drift on the blues density b due to gradients in the reds density r. 
We note that the name "cross-diffusion" to refer to such a term might be confusing, since it 
is really a drift, but this is a common terminology in the literature.— ^ 



C. Self-diffusion 

The self- diffusion coefficient is different to the collective and cross- diffusion coefficients 
in that it is a diffusion coefficient intrinsically attached to an individual tagged particle, 
and which may be related to its mean-square displacement. In contrast, the previous two 
coefficients relate a diffusive flux to the concentration gradient of many particles,— that is, 
to the total concentration of N b or N r particles. Because the self-diffusion is a macroscopic 
property of an individual particle, its analysis in the current framework requires us to tag a 
single particle in the population-level model. We can do this by coloring one particle (the 
tagged particle) in red, N r = 1, leaving the remaining N — 1 particles to be blue particles. 

Setting N r = 1 and N b = N, D b = D r = 1, e b = e r = e br = e, and fj = in ( T22|) . gives 

db 

-(x, t) = V x • [V x 6 +(N- l)t d a bV x b + e d /3 bV x r - e d 7 rV x 6] , (26a) 
— (x,t) = V x • [V x r - iVeVV x r + Ne d f3rVJ)] , (26b) 

where 

2(d-l)7T „ (2d-l)7T 7T , , 
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for d = 2 or 3. Then the self-diffusion coefficient of the tagged red particle is 

D s (b(x, t)) = l- Ne d jb(x, t). (27) 

Hence we find that the self-diffusion coefficient decreases for increasing excluded volume or, 
in other words, that it is reduced relative to point particles. 

Let us compare this to the one-species collective diffusion coefficient. Since the red particle 
is identical to all the blue particles, we untag it so that if the initial densities are the same 
then r = b ( = p, say) and both equations (I26al) and (|26bl) give 

^ (x, t) = V x • ( V x p + Nt d a P V xP ) , (28) 

since /3 — 7 = a. The diffusion coefficient of p is then 

L> c (p(x, *)) = ! + Ne d ap(x, t), (29) 

which coincides with the effective collective diffusion coefficient in (PQ). Thus the collective 
diffusion coefficient D c is increased relative to point particles. This apparent contradiction 
may be understood as follows: the diffusion of any particular particle is impeded by its 
collisions with other particles. However, these collisions bias the random walk towards areas 
of low particle density, so that the overall spread of all particles is faster. When we look 
at the equation for the tagged red ( I26bl) the diffusion is reduced relative to point particles, 
but there is also the drift term due to the gradient in the blues' density. The latter term 
is the dominant one when the blues and reds are the same species, since its coefficient is 
(2d — 1) times larger than the self-diffusion coefficient. The end result is that when the 
two are combined into a single term, the collective diffusion is enhanced relative to point 
particles, as seen in ( f29l) . 

Written in terms of the volume concentration c = (ftp, where = 7rNe d /2d is the volume 
fraction, the self-diffusion coefficient (|27|) and collective diffusion coefficient ( 129]) read 

D s (c) = l-2c } D c (c) = l + 4(d- l)c. (30) 

These expressions are in agreement with the values found in the literature using different 
methods.— 1 ^ Note that the self-diffusion coefficient is independent of the problem dimen- 
sion unlike the collective diffusion coefficient (but note also that it is not defined for one- 
dimensional systems, because the hard-core interaction restricts the allowed motions in one 
dimension^) . 
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D. Experimental measurements of diffusion coefficients 

1. Mean squared displacement (MSD) 

The self-diffusion coefficient of a tagged particle can be described by using the particles' 
mean-square displacement (MSD), defined by MSD(t) = (||Xj(£) — X.;(0)|| 2 ), where X*(t) is 
the position of the zth particle at time t and the angular bracket denotes an ensemble average 
(using that all particles are physically identical). If we keep the convention of coloring in red 
the tagged particle, the following relation is obtained from the second moment of its PDE 
(]26bp (with zero drifts): 

(d/dt) MSD(t) = 2dD s {<j)), 

where d is the problem dimension, D s is given in ( 130]) . and c = <p for the system in equilib- 
rium. This relation, known as the Einstein relation, is more commonly expressed as 

D s ((f)) = lim MSD(t)/(2dt). (31) 

t— >oo 

It relates a macroscopic quantity, the self-diffusion coefficient, with a microscopic quantity, 
the mean-square displacement. The latter can be measured in stochastic simulations of the 
particle system, which we shall use to check the validity of ( 130]) . 

We perform Monte Carlo (MC) simulations of the 2iV-coupled SDE j2| in a two-di- 
mensional periodic box of 400 disks (N = 400) under no external force (f& = f r = 0) and 
uniform initial distribution [Pq(x) = 1]. In order to achieve a better quantitative comparison, 
we employ the event-driven Brownian-dynamics (ED-BD) simulation scheme based on De 
Michele's algorithm. 34 We consider different volume fractions 0, ranging between and 0.1, 
and the particles' size e is chosen to achieve the desired volume fraction. Fig. [2]^a) shows 
an example plot of the mean square displacement of the disks as a function of time when 
= 0.01. Note that MSD(t) varies linearly with t for all times, indicating that the system 
is already in the stationary at the initial simulation time. This is because we have thrown 
away the thermalization transition of our simulations. From its slope at long times the 
self-diffusion coefficient may be extracted using (|3"T]) . Varying the volume fraction in the 
simulation, we obtain points (0, D s ((j))) which are plotted in Fig. Mjo) as red circles. The 
theoretical curve D s ((p), shown as a black dashed-line, compares well with the measured 
values. 
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FIG. 2. Results from stochastic simulation of a two-dimensional periodic system with N = 400 
hard disks and variable volume fraction <p (achieved by changing the particles' diameter e). (a) 
Mean-square displacement MSD(t) for a volume fraction (j) = 0.01. (b) Self-diffusion coefficient 
D s ((j)) for volume fractions of up to 10%. Measured values from stochastic simulations using (|31|) 
(red circles) and theoretical prediction ([30]) (dashed line). 

2. Fluorescence recovery after photobleaching (FRAP) 

Fluorescence recovery after photobleaching (FRAP) is an experimental technique for mea- 
suring the mobility of fluorescent particles. Since the introduction of noninvasive fluorescent 
tagging with fluorescent proteins, this technique has become widely used to study protein 
dynamics in living cells.— In a FRAP experiment, a subregion of the cell (the sample vol- 
ume) is photobleached with a laser beam, causing the molecules contained in it to loose their 
fluorescence irreversibly. As a result, two species are formed, the photobleached molecules 
(inside the sample volume) and the fluorescent molecules (outside the sample volume). Sub- 
sequently, the recovery of the equilibrium from this perturbed initial state is monitored by 
measuring the fluorescence intensity in the sample volume. The resulting curve of intensity 
against time, which can be related to the concentration of the fluorescent species, is then 
used to estimate the overall mobility of the molecule.— 

Several models to fit simulated curves to experimental data have been proposed, most of 
which assume the transport mechanism of proteins to be diffusive.— The standard method 
is based on the work of Axelrod et. alM, in which a linear diffusion equation is used to 
model the fluorescence recovery in a two-dimensional infinite domain. Since the motion of 
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molecules inside the cell can be influenced by many complex interactions such as excluded- 
volume effects or binding events, the estimated diffusion coefficient is often referred to as 
the effective or apparent diffusion coefficient.— 1 ^ 1 ^ Modifications of the diffusion model of 
Axelrod to estimate parameters other than the diffusion coefficient, such as the immobile 
fractions or binding rates, or to account for multiple diffusive species, have also been used.— ^ 
In contrast with the MSD, which is a single-particle tracking method, FRAP is an 
ensemble-averaged method describing the averaged diffusive properties of many fluorescent 
particles.— As a result, even in the simplest scenario of diffusion alone (without binding 
or immobile particles), the diffusion measurements from MSD or FRAP will differ unless 
particles are interaction-free (since one is measuring self-diffusion and the other some for of 
collective diffusion). Nevertheless, FRAP experiments have beed extensively used to mea- 
sure the diffusion coefficient of one molecule, and erroneously seen as an alternative to the 
MSD (for example, when the latter is not feasible due to limitations in labeling or rapid 
diffusion) .— 

Therefore, care should be taken to interpret the diffusion coefficient estimated from 
FRAP experiments. The first point to note is that the FRAP diffusion coefficient can- 
not be identified as the self-diffusion coefficient of the molecule, often termed the anomalous 
sub diffusion. 38 Instead, it characterizes the mobility of the whole fluorescent population. 
Secondly, the standard linear diffusion equation ignores interactions between the photo- 
bleached and fluorescent species. Is not clear whether a "pure" diffusion measure can be 
obtained from the FRAP experiments or, as we have seen in Sec. IIII B\ it is a combination of 
collective diffusion and drift due to gradients of the photobleached species (cross-diffusion) 
that it is in fact measured. 

The two-species model ( 122]) can be used to model the FRAP experiment. For example, 
consider the simple setting a two-dimensional unit square domain, pure diffusion and a circu- 
lar photobleached area (sample area) of radius w,— with e<w<l. The system is initially 
in equilibrium, so that the concentration is uniform in Q and equal to the volume fraction 
= 7cNe d /2d. A laser of intensity / causes the portion of particles contained in ||x|| < w 
to be photobleached, which we identify as the blue particles. The fluorescent particles com- 
prise the red species. We consider volume concentrations rather than probability densities 
to relate them to FRAP measurements. The number of bleached particles is Nf, = irw 2 N 
and their initial volume concentration is c&(x, 0) = <ft for ||x|| < w and otherwise. For the 
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fluorescent particles, N r = N(l — tcw 2 ) and c r (x, 0) = (ft f° r Il x ll > w - Since initially all 
particles are identical, they have equal diffusion coefficient D. Using that N r is large such 
that N r — 1 « N r , the equation for the fluorescent (red) species (122bj) reads 
Oc 

-^(x, t) = DV X . {[1 + 4(d - l)c r - 2c b ]V x c r + 2(2rf - l)c r V x c fe } , (32) 

where c r = 7rN r e d r /2d. An analogous equation is obtained for the bleached (blue) species. 
Solving this system of equations for c r and c b with the initial conditions described above until 
equilibrium is reached [uniform concentrations c r (x, oo) = (1 — nw 2 )(ft and Q,(x, oo) = nw 2 (ft}, 
the theoretical recovery curve, which is related to the integral over the sample area of 
J(x)c r (x, t), 35 could be compared with the experimental recovery curve. 



IV. BASIC PROPERTIES 

In this section we discuss some basic properties of the cross-diffusion model (122 p . such 
as its free-energy functional, equilibria and stability of solutions. We restrict ourselves to 
the case when the force fields are potential forces, that is, f&(x) = — V x Vb(x) and f r (x) = 
— V x VJ.(x); i.e., we consider the system 

— (x, t) = V x • (A, [1 + (N b - l)e d b ab] V x 6 + V x H(x)& 

dt V (33a) 

+ N r e d br {D b % 6V x r - lb rVJ>] + 7fe V x [K( x ) - H( x )] &r}) , 

— (x,t) = V x ■ (D r [I + (N r - l)e d r ar] V x r + V x K( x )r 

dt V (33b) 

+ N b e d br {D r \p r rV x b - lr bV x r] + 7r V x [H( x ) - K( x )] &r}) , 

in Q x ]R + with no-flux boundary conditions on dfl x M + and initial conditions 

&(x,0) = &o(x), r(x,0)=r (x) in fi. (33c) 
The coefficients a, and 7i (z = 6, r) are all non- negative and given in (122dl) . 



A. Gradient flow structure and free energy 



In Ref. 
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we found that, when f(x) = — V x l / (x), the equation for one species (pP) can 
be written as gradient flow^ 



dp 
dt 



+ V x • (pu) = 0, 



(34a) 
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with u = — V X [.D \ogp + 2Dctd(N — l)e d p + V(x)]. (In Ref. [ill u an d -^(p) were defined as 
— u and -F(p) respectively.) Here u can be thought of as a "flow" down the gradient of the 
free energy J-(p) associated to equation flTJ, u = — V x ^p, with 



J r {p)= / D \p\ogp + a(N — l)e d p 2 ] dx + / K(x)pdx. 



(34b) 



The first integral corresponds to the internal energy, which increases with excluded-volume 
effects, and the second integral is the potential energy. Then the free energy is non increasing 
in time when evaluated along a solution of ([T]). 4 - The gradient flow structure fl34l) is very 
useful since it brings more tools to study the trend to equilibrium,— with the free energy 
functional ( I34bj) "encoding" all the properties of the flow. 

The scalar gradient flow structure f)34ap becomes, in the case of two species, 12 




MV, 




(35) 



where J 7 = J-"(6, r) is, again, a scalar free-energy functional and M = M(6, r) is a two- 
dimensional matrix denoted the mobility matrix (which must be positive semi-definite from 
the definition of free energy). In this section we examine under which conditions the two- 
species system (133]) admits an explicit gradient flow representation (135]) valid to 0(ef,ef), 
bearing in mind that eb,e r ,eb r <C 1. It should be pointed out that the transformation of 
( I33ap - (133bj) into the structure (135]) is not straightforward in general. In the following we 
present two situations in which our system admits a gradient flow structure: for a large 
number of particles (Subsec. IIV A II) and when the drift terms become zero (Subsec. IIV A 21) . 



1. Large number of particles approximation 

Assume that the number of blue and red particles is large such that A^ — 1 ps Aj, and 
A" r — 1 w N r , and that the two species have the same diffusion coefficient (which we can set 
to unity without loss of generality). Then the system ( 133]) for the densities b and r can be 
rewritten in terms of the number densities b = N^b and f = N r r in the gradient flow form 
( 135]) . with free-energy functional 

JF(b, r)= J b log 6 + f log f + bV b + rV r + | (e d b 2 + 2t d hr br + t d f 2 ) dx, (36a) 
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and mobility matrix 



Mft^l^f' / 4 ,.J, (36b) 



where a = 2(d — l)it/d and 7 = tt/g? is the simplified version of 7, (i = b,r) in (122dj) when 



the diffusivities of blues and reds are equal. [The coefficient /3j disappears in Eq. f l3T)|) using 
that /3j = a + 7 when D b — D r — 1.] The mobility matrix M is positive definite if and only 
if 

T e£(& + r)<l, (37) 

using that b,f > 0. This upper bound on the total number density b + f gives a limit 
of validity of the model, and must be satisfied pointwise in Q. Nevertheless, to get an 
approximate idea of the upper bound on volume fraction 0, assume that the concentrations 
are constant, b = N b and f = N r and that e b = e r . Then, using that the volume fraction 
is = ld e trX^b + N r ), Eq. ( J37j) becomes 4> < 0.5. Therefore, we see that our model for 
two species in the case of large number of particles and equal diffusivities would break down 
when the volume fraction reaches one half. 



2. Zero potential 



Now consider the system (133]) with zero potential, Vj, = V r = 0. In this case our problem 
has also a gradient flow structure (!35l) with 



T(b,r) 



Mogfc + rlogr + ^{N b - l)e d b ab 2 + ^{N r - l)e d r r 2 + 6 e d br br 

_ 



dx, 



(38a) 



and 



M(b,r) 



f D b b{l - N r e d hrlb r) D b e d br (N r (3 b - 6)6r^ 
K D r ei{N b p r - Q)br D r r{l - N b e d br lr b) y 



(38b) 



where B is a free parameter. There are two relevant cases for the gradient flow structure 
( |38l) depending on the value given to 6. First, the free energy (I38ap can be chosen to be 
the sum of the two one-species entropies like ( I34bl) by setting = 0. Second, the mobility 
tensor can be adjusted to be symmetric setting 9 = (N b D r f3 r — N r D b (3 b ) / \D r — D b ) provided 
D b ^ D r . (In case D b = D r , we can still obtain a symmetric mobility matrix while at the 
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same time setting = by rewriting (138]) in terms of the number densities b and f as in 
the previous subsection.) The determinant of the mobility matrix (138bl) is 

det Mft r) = D„DM (l - f ^g^ *) + 

where note that the free parameter G only enters at O(e^). The mobility matrix is positive 
definite if (neglecting the C(e^) terms) 

2tt ( A-Nfc& + D b N r r) 



d D b + D T 



< r < 1, (39) 



using that b, r > 0. As before, if we assume that b = r = 1 (we have reached the equilibrium), 
D b = D r and e b = e r , the upper bound f[3"9"j) becomes < 0.5 as in the previous subsection. 
The point at which the mobility matrix becomes singular can be related with the stability 
of the equilibrium states of the system, see Sec. IIV CI 



B. Equilibria 

We compute the stationary solutions of the nonlinear diffusion model (155|) . which we 
denote by b s and r s . Note that, in the case of point particles (e& = e r = 0), the equilibria 
are trivial as the system for b s and r s decouples and we find 

6,(x) = C b exp[-V b (x)/D b ], r s (x) = C r exp[-V r (x)/D r ], (40) 

where C b and C r are the normalization constants, i.e., Ci = J n exp[— Vi{x)/Dj\ dx. For 
finite-size particles, we first consider the cases for which we have found a gradient flow 
structure (1351) in Sec. IIV Al and then we consider the general case in two ways: solving the 
steady states of (|33|) with a finite difference numerical scheme, and similarly the stationary 
distribution of the particle description (j2J) with stochastic simulations using the Metropolis- 
Hastings (MH) algorithm.— 



1. From the free-energy description 



One big advantage of the system with a gradient flow structure is that the steady states 
can be found in a straightforward manner by minimizing the free energy J 7 . If the mobility 
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matrix M(6 s ,r s ) is positive definite, then finding the steady states of (155]) is equivalent to 
finding constants Xb and Xr such that 

d h JF{b s ,r s ) = Xb, d r JF(b s ,r s ) = Xr, (41) 

with J b s = J r s = 1. Note that the no-flux boundary conditions on dfl are automatically 
satisfied by imposing d b T and d r T to be constant. 

In the case of no external potentials (Sec. IIV A2p . we see that the obvious constant states 
b s = r s = 1 are the only stationary solutions of (|38|) . In the case of a large number of 
particles (Sec. IIV A lj) . the stationary solutions b s and f s (recall we used number densities 
for this case) are found by imposing (|4T!) on the associated free energy (I36al) . Then the 
problem reduces to 

\ogb + V b + a(efb + e d br r) = X b, , s 

( 42 ) 

\ogf + V r + a(efr + eib) = X r, 
such that f n b s = N b and J n r s = N r . This problem can be solved with the Newton-Raphson 
method, discretizing Q into J grid points, approximating the normalization integrals with a 
quadrature, and solving for 2 J + 2 unknowns (xb, Xr, b\, and r k , for z = 1, . . . , J). 



2. From the particle system via the Metropolis-Hastings algorithm 

As we did with the time-dependent case, we now compare the stationary solution of the 
reduced macroscopic model (l33l) to the stationary distribution of the full microscopic model 
( )4al) . Under external potential forces, ([3]) becomes, on rearranging, 

dp - 

— (x,t) = AV Xi • {P V Xi [logP + ViteVA]} , (43) 

where = _D b , = Vf, for z < Afc and A = A, V% = K otherwise. Denoting D the 

A^— diagonal matrix with diagonal entries A, we can write 

dP 

— + VT©(P«) = 0, (44) 

where u = — V^flogP + J2iLi ^i( x i)/A] is the "flow velocity" vector field. Using that D is 
non-singular and that no-flux boundary conditions hold, the stationary solution P s of ( )44l) 
is P s (x) = Cexp[— T-L(x)), where %{x) is given by: 

, y?, Vi(xi)/Di, xefl", 
U{x) = { ^ l=1 K }l ' e ' (45) 

oo, otherwise. 
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3. Examples 



Two examples of stationary densities are shown in Figures [3] and [5j The stationary 
solution of Eq. fl33|) . computed using (J12|) . is compared with stochastic simulations of the 
stationary distribution of the full iV-particle SDE (J2J) with the MH algorithm. 




-0 5 ^ 5 -0 5 ^ 0.5 -0.5 0.5 -0.5 0.5 

FIG. 3. Stationary marginal densities b s (x) ([I— IV]-b) and r s (x) ([I-IV]-r) for point and finite-size 
particles, with V b = Wy, V r = 5y, D b = D r = 1, N b = 600 and N r = 200. (I) Solutions b s (x) and 
r s (x) of (Ti2|) for point particles (e b = e r = 0). (II) Histograms for e b = e r = 0. (Ill) Solutions b s (x) 
and r s (x) of (|42p for finite-size particles (e b = 0.01, e r = 0.015). (IV) Histograms for e b = 0.01, 
e r = 0.015. Histograms computed from 10 steps of the MH algorithm. All four plots on the left 
and right have respectively the same color bar (note inverted color scale for reds and blues). 



Fig. E] corresponds to the stationary solution under a 'gravitational' force in the direction 
-e y , with n = [-1/2, -1/2], e b = 0.01, e r = 0.015, D b = D r = 1, N b = 600, and iV r = 200. 
We suppose that the blue particles (four plots on the left) are heavier than the red particles 
(four plots on the right), and that therefore they feel a stronger force, i b = — 10e y versus 
f r = — 5e r While both blue and red particles accumulate at the bottom when finite-size 
effects are ignored [Fig. Ell-b) and Fig. E](I-r)], the blue particles accumulate at the bottom 
[Fig.[3](III-b)] and force the red particles upwards [Fig. [3](III-r)] when they are not (note there 
is zero probability of finding a red particle at y = —0.5). The averages of b s and r s across x 
are shown in Fig. HI The agreement between the model (|42p and the stochastic simulations 
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is good in all cases, except in the region near y = —0.5 for the red finite-size particles 
[compare Fig. |3tlII-r) and Fig. EflV-r)]. A possible explanation for this disagreement is that 
the variability of r s near that boundary occurs in a region of size equal to the size of the 
histogram bins. 
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FIG. 4. Averaged stationary densities across x, {b s ) x and (r s ) x , corresponding to the 8 cases shown 
in Fig. [3l Comparison between stationary solutions of (|42p (curves) and histograms obtained from 
MH simulations (data points). 

Fig. corresponds to the stationary solution under a symmetric bivariate Gaussian po- 
tential of the form 

£(x; fx, a 2 ) = 1/(2trx 2 ) exp{-[(x - /i) 2 + (y - ^) 2 }/2a 2 }. 

The parameters are V b = -0.1£(x; 0, 0.05), V r = 0.5£(x; 0.2, 0.05), e b = e r = 0.02, D b = 
D r — 1, N b — N r — 400 and Q = [—1/2,1/2]. For point particles, the stationary solu- 
tions preserve the radial shape and centers of their respective potentials Vb and V r , i.e., 
b s oc e~ Vb and r s oc e~ Vr [Fig. [5]^I-b) and Fig. EJl-r)]. However, we can appreciate the dis- 
torted/asymmetric shape of the reds' density r s when size effects are included [Fig. E^III-r)]. 
Also, in the blues' density we can observe clearly how there is a competition between the 
potential well and the finite-size repulsion — the particle density b s inside the well is reduced 
for finite-size particles. Again, the agreement between the model (T42|) and the stochastic 
simulations is excellent. 
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FIG. 5. Stationary marginal densities b s (x.) ([I-IV]-b) and r s (x) ([I-IV]-r) for point and finite-size 
particles, with V b = -0.1£(x; 0, 0.05), V r = 0.5G(x; 0.2, 0.05), D b = D r = 1 and N b = N r = 400. 
(I) Solutions of (|4"2|) for point particles (e;, = e r = 0), b s oc e~ Vb and r s oc e~ Vr . (II) Histograms for 
e b = e r = 0. (Ill) Solutions of (j4~2j) 6 s (x) and r s (x) for finite-size particles (e;, = e r = 0.02). (IV) 
Histograms for e b = e r = 0.02. Histograms computed from 10 7 steps of the MH algorithm. All four 
plots on the left and right have respectively the same color bar. 

C. Linear stability 

It is of interest to compare the upper bounds obtained from the gradient flow structure in 
Sec. IIV Al (looking when the mobility matrix M becomes negative definite) with a classical 
linear stability analysis. We consider a simple case here but the analysis is straightforward 
for more general cases. 

Consider the system ( 133|) with linear potential forces, that is, potential forces of the form 
Vf,(x) = Vfi-x and V r (x) = v r -x. In such cases the equilibrium states are simply b s = r s = 1. 
We make the following linearization around the equilibrium states, 

b = 1 + 5A b exp (at + zk • x) , r = 1 + 5A r exp (at + ik • x) . 

Inserting these into and neglecting 0(8 2 ) terms yields a system B(cr, k) ) =0. The 
condition for a non-zero solution, det B = 0, is the dispersion relation. For the basic case 
e b = e r = e, D b = D r = 1 and v& = v r = 0, we find that one solution of det B = is always 
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negative and the other one is 

<r(k) = ||k|| 2 |-1 + tj [2{d - 1) + N b + N r ] | . (46) 

This corresponds to a zero- wavelength infinite growth rate instability when <p+e d 7r(d— l)/d > 
1/2, where <fi is the particle volume fraction. 

The condition that o < in f|46|) (for stability) is equivalent to the condition fl37|) in 
Sec. II V ATI found from the mobility matrix under the assumption Nj,, N r ^> 1. Observe that 
in the condition f )37|) the magnitude of the drift terms did not play a role. The numerical 
exploration of sgn(<r) for several drifts confirms this: an arbitrarily large drift cannot change 
the sign of o. We emphasize again that this instability represents a breakdown of the model 
reduction, not a true instability in the original particle system. 



D. Symmetrizability of the system and the Onsager relations 

In Sec. II III we have seen how our multicomponent diffusion system involves the study of 
a diffusion matrix that describes how the flux of one component is influenced by its own 
density gradient and the density gradient of the other component in the system. These ideas 
can be related with the thermodynamic Onsager reciprocal relations, which establish that in 
a system fluctuating around its equilibrium a relationship between certain fluxes and ther- 
modynamic forces must hold. 43 The Onsager relations are defined assuming the fluctuations 
around the equilibrium are small (so that the relationship is linear). In particular, if in a 
system we have the following relations between fluxes Jj and forces X i: 

J l = Y,L tk X k , (47) 

k 

the Onsager reciprocal relations requires symmetry in the cross-terms, that is, = L^i- 
The Onsager relations are a macroscopic consequence of microscopic time reversibility (prin- 
ciple of detailed balance).™^' 35 - ) Note that coefficients can be nonlinear functions of the 
variables.— ( p ' 64 ^ Gupta and Cooper— study the relationship between the matrix L and the 
diffusion matrix D in a linear multicomponent diffusion, and find that D must be positive 
definite for the Onsager relations to hold. It should be pointed out that while the Onsager 
relations have been named the fourth law in thermodynamics by some authors, their validity 
has yet to be indisputably established. For instance, many valid multicomponent diffusion 
models have been found to be inconsistent with these relations. 30 
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We proceed to investigate in which cases, if any, our cross- diffusion model (1231) is consis- 
tent with the Onsager relations. It is appropriate to consider the free-energy gradient with 
respect to the densities to be the force X driving the system to its minimum free-energy 
equilibrium stated 13 ' 35 ^ The gradient-flow structure (I31)|) fits with the form required in (|47[) . 
and the question is for which cases the mobility matrix M(6, r), L in (j47p . is symmetric. In 
Sec. IIV Al we found two situations for which our system satisfies this: the mobility matrices 
for the large number of particles' approximation f!36bj) and for a zero potential f!38bj) (with 
the appropriate choice of the parameter 0) are indeed symmetric. 

The fact that in a system with a positive definite diffusion matrix the Onsager relations 
hold may be related to analytical work on parabolic systems. It is well known that, in 
hyperbolic or parabolic systems, the existence of a free-energy functional is equivalent to 
the existence of a change of unknowns which "symmetrizes" the system.—^ For parabolic 
systems, "symmetrization" means that the transformed diffusion matrix is symmetric and 
positive definite.^ Note that in Sec. IIV Al the change of unknowns consisted of going from 
(b, r) to the so-called entropy variables (df,E, c^-E 1 ).— The equivalence between Onsager rela- 
tions and symmetrization comes from the fact that the symmetry is a necessary condition for 
the entropy production rate to have a sign, itself a condition for the system to be compatible 
with the second law of thermodynamics. Although these analytical results seem promising 
in order to find a free energy for the general form of our system with non-zero potentials 
( l33l) . it should be pointed out that finding the change of variables that make our system 



symmetric (in the sense described in Ref. |47J) can be in general as challenging as finding the 
free energy itself. A first step would be to find a change of unknowns for which the system 
can be put in a form with no drift terms. 
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that the original diffusion 



To conclude this section, we check that the result in Ref. 
matrix must be positive definite holds in our case. To 0(ef r ), our diffusion matrix f !23bp has 
eigenvalues 

Xi = D { + ^D^d-im - l)e?i(x) - ^-^VtiW) , (48) 

for i = b, j = r, and vice versa. A lower bound is A« > Di + ^jD^Niefi — Nje'jj}, from which 
we find that A r > (since we must have small volume fraction, i.e., Nbef + N r ef < 1). 
Therefore, provided there is a small volume fraction, our diffusion matrix is positive definite 
and hence, as we have already mentioned above, the Onsager relations hold for our system 
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with zero-potentials. 



V. CONCLUSIONS 

In this paper we have considered the diffusion of two interacting species of hard spheres, 



extending the model derived in Ref. [ll| to incorporate particles of different sizes, different 
diffusivities and under different external forces. The result is a nonlinear cross-diffusion 
PDE system for the two marginal probability densities associated to each species. This 
approach enables us to describe a complicated system of interacting particles with a simple 
the continuum PDE model whilst capturing the hard-core interactions at the particle level. 
These interactions emerge as nontrivial nonlinear terms in the continuum model, involving 
cross-coupling terms which can be interpreted in terms of the inter-species competition at 
the population- level. In addition to providing some insight on the system's behavior, the 
continuum model is relatively easy to solve and analyze. We have assessed the validity of 
our continuum approach to predict the behavior of the system by comparing its numerical 
solutions with stochastic simulations of the discrete particle-based model. We have found 
very good agreement between the two, supporting the idea that by solving a simple system 
of PDEs we can capture the same population-level behavior observed after many repetitions 
of expensive stochastic simulations. 

Our analysis is valid in the limit of small but finite particle volume fractions so that 
pairwise interactions are the dominant ones. This excludes situations close to the jamming 
limit. Our method uses matched asymptotic expansions in the particle volume fraction to 
derive the continuum model in a systematic way as a perturbation of the interaction-free 
case. Because of the perturbative nature of the method, we expect its accuracy to decline 
as the volume fraction increases and eventually cease to be valid. In particular, by writing 
the system in a gradient-flow form in terms of the free energy functional and a mobility 
matrix, we can use the singularity of the mobility matrix as an indicator of the break-down 
of the model; this occurs at roughly 50% volume fraction. Despite the limitation of a low- 
volume fraction, we believe our method can provide insight into the mechanisms by which 
particle-level characteristics emerge at the population-level. 

Our method is not based in the thermodynamic limit which requires the number of 
particles iV to tend to infinity together with the system volume.— Therefore, the continuum 
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Fokker-Planck (FP) model derived here should not be misinterpreted as a deterministic 
model for the concentrations of the two species in the system, valid only when the number 
of particles N is large. While it can be used in this situation if required, we emphasize that 
in our work the reduced FP model is valid for any iV (one could set N^, N r = 1). In other 
words, the continuum system is not a deterministic model, but rather a PDE system for the 
probabilities of finding a blue and a red particle at a given position and given time. 

Our two-component drift-diffusion model captures an enhancement of the collective dif- 
fusion rate due to excluded-volume interactions between particles of the same species (as 
we had already found in our previous workJA), as well as a reduction of the collective dif- 
fusion due to interactions with particles of the other species. This structure is useful not 
only to study the collective diffusion in terms of the particles' volume fraction, but also to 
analyze the self-diffusion coefficient. The latter describes the evolution of a single tagged 
particle, and it can be extracted from the model by choosing one of the species to have only 
one tagged individual. In contrast to the collective diffusion, which increases with volume 
fraction (or relative to point particles), the self-diffusion coefficient decreases with volume 
fraction. Thus the two species model can be used to characterize transport properties of 
a system of identical particles by distributing them in two subpopulations of iV — 1 and 
one particles respectively. To our knowledge, such a continuum model capable of explaining 
both the collective and individual diffusion coefficient has not been reported in the liter- 
ature/previous works. We have briefly described two experimental procedures to obtain 
diffusion measurements from real systems, namely a single-particle tracking method to mea- 
sure the particle's mean-square displacement (MSD) and the ensemble-averaged method 
FRAP (fluorescence recovery after photobleaching) . While it is well understood that the 
MSD can be related to the self-diffusion coefficient, there is some debate and confusion 
over the interpretation of FRAP experiments.— Several methods based on fitting curves to 
FRAP experimental data exist,— but none appears to have identified the fact that FRAP is 
measuring a collective transport property. We believe that our two species model could be 
used to model the FRAP experiment systematically and provide an improved framework to 
interpret its results. As can be seen from the FRAP model ( 13"2|) . in general, the measured 
quantity will not be a pure diffusion coefficient (such as the collective diffusion coefficient, 
as one might be tempted to think from the fact that FRAP gives an ensemble-averaged 
measurement) but a mixture of collective diffusion with advection due to gradients in the 
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concentration of the photobleached species. 

We have investigated for which values of the parameters the cross-diffusion system accepts 
a gradient-flow structure in terms of a free-energy functional; this structure is useful to study 
the equilibria of the system as well as its stability. Namely, the stationary solutions of the 
system correspond to the minimizers of the free energy, and the stability can be studied 
from the convexity of the free energy functional near its equilibria. 

Previous work on the diffusion of two species with size-exclusion interactions using a 
lattice-based modet^i 8 - led to a continuum population-level description which is different 
from our reduced model (which does not restrict the motion of particles to a lattice). In 
other words, the two approaches (on- and off-lattice models at the particle-level) result in 
different reduced models, even though they are trying to describe the same problem. It 
would be interesting to study which rules one needs to prescribe on the lattice model in 
order to achieve a certain population-level description. We will address this issue in future 
work. 

The model presented in this work can be extended to consider the diffusion of finite-size 
particles through obstacles, which has many important applications in porous media and 
diffusion in biological systems. This may be achieved by setting the diffusivity of one of the 
species (the obstacles) to zero. An advantage of this approach is that it makes it very easy to 
study diffusion through spatially varying concentrations of obstacles. 48 Another interesting 
extension is to consider anisotropic particles 4 - and examine how the continuous PDE model 
changes with nonspherical particles. 
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